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Abstract 


A computationally efficient 3DHZETRN code capable of simulating High charge (Z) and Energy 
( HZE ) and light ions ( including neutrons ) under space-like boundary conditions with enhanced 
neutron and light ion propagation was recently developed for simple shielded objects. Monte 
Carlo (MC) benchmarks were used to verify the 3DHZETRN methodology in slab and spherical 
geometry, and it was shown that 3DHZETRN agrees with MC codes to the degree that various MC 
codes agree among themselves. One limitation in the verification process is that all of the codes 
(3DHZETRN and three MC codes ) utilize different nuclear models/databases. In the present 
report, the new algorithm, with well-defined convergence criteria, is used to quantify the neutron 
leakage from simple geometries to provide means of verifying 3D effects and to provide guidance 
for further code development. 

Introduction 

Early nucleon transport code development relied on Monte Carlo (MC) methods using the straight-ahead 
approximation [Wright et al 1969]. Further development of a full 3D MC code with the Bertini intranuclear cascade 
model was used to test the straight-ahead approximation for which resulting dose and dose equivalent were within 
the statistical uncertainty of the full 3D MC code [Alsmiller et al. 1965]. Such 3D MC codes on computers of that 
time required vast resources even for small sample sizes [Wright 1966, Alsmiller 1967] and encouraged the 
development of more efficient means of radiation exposure evaluation. A 3D MC code that sampled a Bertini model 
derived database giving it much greater computer efficiency was extended to high energies at the NASA Langley 
Research Center (LaRC) [Lambiotte et al. 1971] for atmospheric and space radiation studies but still required 
extensive computer resources for realistic calculations [Foelshe et al. 1974]. As a result, NASA LaRC began a 
search for methods of direct solution of the Boltzmann equation. It was found that proton transport based on direct 
solution of the straight-ahead Boltzmann equation resulted in an efficient code able to reproduce the straight-ahead 
MC using the same database [Wilson and Lamkin 1974, 1975; Larrikin 1974; Wright et al. 1969]. As a result, NASA 
LaRC began a long-term investigation of direct solution of the Boltzmann equation as an alternative to MC methods 
[Wilson and Lamkin 1974, 1975]. The 3D MC codes would play an indispensable role in these developments as a 
test for this alternate evaluation scheme. 

A major element in developing Boltzmann solution methods is the examination of adequate computer 
methods, as was first taken up in the straight-ahead approximation following the MC studies of Wright et al. [1969] 
and Alsmiller et al. [1965]. Both numerical convergence studies [Lamkin 1974] and formal convergence studies 
following Wilson and Lamkin [1974, 1975] were commissioned. The next formal convergence study examined the 
straight-ahead approximation itself as a first order asymptotic expansion in terms of angular factors for which it was 
shown that errors were on the order of the square of the ratio of the beam divergence (few cm) to the radius of 
curvature of the object (few meters in human rated space systems) [Wilson and Khandelwal 1974]. Additional 
analysis using the Fredholm version of the Boltzmann equation, along with asymptotic expansion in angular width, 
demonstrated that the straight-ahead approximation is more accurate than assuming that the High charge (Z) and 
Energy (HZE) fragments proceed with the same velocity as the projectile ion [Wilson 1977]. Neglect of these higher 
order terms resulted in reduction of the Fredholm version to a Volterra equation that could be solved using marching 
procedures [Wilson and Badavi 1986]. 

Solving the Volterra equation as a perturbation series (Neumann series) allowed direct comparison to 20 Ne 
attenuation experimental data obtained at the Lawrence Berkeley Laboratory's Bevatron [Wilson et al. 1984] 
revealing the need for an improved nuclear fragmentation model. The initiation of the NUCFRG model allowed 
improved comparison with atmospheric air shower data [Wilson et al. 1987a,b] and with the 20 Ne Bevatron data 
[Shavers et al. 1990, 1993]. With this success in transporting HZE ions, only the treatment of light ions and neutrons 
required advancement as a formalism compatible with the HZE marching procedures. 

Unlike the HZE fragments that are produced at velocities nearly equal to the projectile ion, neutron and 
light ion induced reactions produce a broad spectrum of secondary particles and cannot be adequately simplified. 
This complication was addressed by Wilson et al. [1988a,b] using perturbative methods wherein a closed form 
solution is first obtained over a small spatial step h and repeated to a depth x = Nh. The approach provided a 
computationally efficient solution over the whole spatial and energy domain. Such marching procedures are fully 
compatible with the HZE code and the integration of the baryon transport code (BRYNTRN) into the HZE code 
provided the initial HZETRN version [Wilson et al. 1991]. Cucinotta [1993] added the remaining light ions 
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(deuteron, triton, helion, and alpha) to BRYNTRN and integrated the revised code into HZETRN. Code verification 
with analytic solution methods were accomplished by Wilson et al. [1987b], Ganapol et al. [1991], Lamkin [1994], 
and Shinn et al. [1998]. Code validation was summarized by Shinn et al. [1998]. Later work completed by Slaba et 
al. [2010a] greatly improved overall code efficiency and reduced numerical discretization error associated with the 
marching algorithms. Blattnig et al. [2004] and Norman et al. [2012, 2013] also integrated a coupled model for pion 
and muon transport and the associated electromagnetic cascade. These efforts largely completed the development of 
HZETRN in the straight-ahead approximation. Further evaluation of the accuracy of the straight-ahead formalism 
requires formal improvements with added detail to the description that can be taken in two directions: corrections to 
the straight-ahead approximation or addition to the physical process descriptions within the HZETRN formalism. 

An obvious first correction relates to the assumption that even particles produced in the backward 
hemisphere travel forward in the straight-ahead approximation. Hence as a correction, the backward hemispheric 
produced neutrons were assumed to travel straight backward as the next level of development (bi-directional 
approximation). Note, this changes the character of the solution as forward and backward propagation is coupled in 
the nature of the Fredholm equation. Clowdsley et al. [2000, 2001] solved the forward/backward neutron equations 
in the multigroup approximation and using collocation methods. Heinbockel et al. [2003] compared additional 
solution methods among themselves. Slaba et al. [2010b] fully coupled the forward and backward neutron transport 
equations and improved numerical efficiency by casting the system in a linear algebraic equation format (as a 
discrete version of the Fredholm equation). Extensive comparisons have been made with available MC codes 
[Heinbockel et al. 201 la, b; Slaba et al. 2011, 2013] and found to agree with the various MC codes as well as the MC 
codes agree among themselves. 

The next logical step allowing an objective evaluation of the bi-directional approximation error is to allow 
fuller angular dependence within the Fredholm formalism. A first step was implemented in a perturbative approach 
in finite objects and compared to MC results [Wilson et al. 2014a, b,c, 2015]. This evaluation of 3D effects by 
comparing with MC evaluations met with some of the limitations of the various MC codes which use differing 
nuclear models and databases resulting in large disparity [Heinbockel et al. 201 1 a, b] so that 3D effects are 
somewhat obscured [Wilson et al. 2014a, b]. The nuclear cross sections enter in two ways; first, they differ in the 
overall production of neutrons/light ions among the transport models, and second, they differ in attenuation and 
scattering cross sections as particles move through materials. Herein, means of reducing the effects of the differing 
neutron production terms among various transport models will be used to more unambiguously evaluate 3D effects 
in neutron transport. 

Deterministic Code Development 

The following two sections summarize deterministic code formalism and developments resulting in the 
HZETRN and 3DHZETRN codes. Much of the information has been documented elsewhere (see in particular 
Wilson et al. [2014a, b,c, 2015] and references contained therein) but is being included here for completeness. 

The relevant transport equations are the linear Boltzmann equations derived on the basis of conservation 
principles [Wilson et al. 1991, 2005] for the flux (or fluence) density <p j (x : 12, E) of a j type particle in the 

continuous slowing down approximation (CSDA) in which atomic processes are described by the stopping power 
Sj(E) for each ion type j (vanishes for neutrons j = n) as 


B [^(xME)} = ZJ7L a^{E,E\n,n')4> k (x,n',E')d£l'dE' , (1) 

k 

where the differential operator on the left hand side is defined as 

B[4>.(x,n,E)] = Sl.V4>.(x,n,E) - j--E[ Sj {E)4>.(x,n,E)] + a^E^xME) . ( 2 ) 

Equation (1) is solved subject to a boundary condition over the enclosure of the solution domain as shown in Fig. 1. 
The terms cr. (E) and <j j k (E,E\EL,EL } ) are the media macroscopic cross sections and include nuclear elastic and 
reaction contributions and have been discussed in more detail elsewhere [Wilson et al. 1991, 2005]. 
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Fig. 1. Geometric relations of quantities useful in solving equation (1). The symbol n represents a unit normal 

vector. 


Application of CSDA in both laboratory and space shielding has been wide and the resulting errors 
discussed elsewhere [Wilson et al. 1984, Shavers et al. 1993, Tweed et al. 2006a, b]. Equation (1) can be rewritten as 
a Fredholm equation [Wilson 1977] given by 


^,n, £ )- S ;W 


Sj(E)Pj(E) 


[r(n,*),n,E 


+- 
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E f E E /“ f 4 w p j > jk (E',E",n,n ')c/> k (x n ,si\e ")< m 'dE"dE\ 


( 3 ) 


SAE)PAE)^Je 


with 


£ 7 = Rpp-d + R^E)}, (4) 

x sl =x + [R j (E)-R j (E')]n. (5) 

In equations (1) and (3), r(Q,jc) is the point on the boundary connected to x along -O; p is the projection of x onto 
O (Fig. 1); d is the projection of r(Q,jc) onto O; Rj(E ) is the average distance a type j ion with energy E will travel 
before losing all of its energy to excitation/ionization of atomic electrons, and Pj(E) is the probability a type j ion 
with energy E will not have a nuclear interaction in coming to rest in the media. The range-energy relation is given 
by 


R 3 (E) 


jl . 

jJ o S j (E ') 


and the nuclear attenuation function is given by 


Pj(E) = exp 


A p zjEl! 

o Sj(E') 


dE } , 


( 6 ) 


( 7 ) 


where o] includes nuclear elastic and reactive processes. The zero charge limit of equations (1) and (3) is discussed 
elsewhere [Wilson et al. 1991]. One obstacle to solving either equation (1) or (3) is the need to evaluate the integral 
over d£l' at arbitrary locations within the media and development of methods to handle this limitation expediently in 
ways amenable to efficient computational practice. 

Atomic interactions limit the contributions of charged particles in the transport process. For example, the 
protons and alpha particles produced in aluminum below 100 MeV/n contribute to the fluence only within a few 
centimeters of their collisional source and the heavier ions are even more restricted. This is an important factor in 
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that the transported secondary charged particle flux tends to be small at low energies and the role of additional 
nuclear reactions by these low energy ions are likewise limited, an important fact to be implemented later in the 
current development. 

As noted above, a prime limitation in the evaluation of equation (1) or (3) is evaluation of the integral over 
d£l' at any arbitrary location within the media. The approach to a practical solution of equation (1) or (3) is to 
develop a progression of solutions from the simple to increasingly complex allowing early implementation of high- 
performance computational procedures and establishing a converging sequence of approximations with established 
accuracy criteria and means of verification and validation. The first step leading to the lowest order solution reduces 
the evaluation by introducing the straight-ahead approximation as guided by the nucleon transport studies of 
Alsmiller et al. [1965] using Monte Carlo methods in which the differential cross sections were approximated as 

a jk {E,E',Sl,Sl') = a jk (E, E ')«(« -O'). (8) 

resulting in dose and dose equivalent per unit fluence to be within the statistical uncertainty of the Monte Carlo 
result obtained using the fully angle dependent cross sections in slab geometry [Wilson et al. 1991]. 

In the early Boltzmann transport studies of Wilson and Lamkin [1974, 1975] and Lamkin [1974], it was 
demonstrated that after neutrons are produced by a proton beam, the recoupling of neutrons back to the proton field 
is limited, but the neutron fields tend to build-up with increasing penetration depth. This mainly occurs since many 
protons produced by neutron-induced reactions are of lower energy and are quickly removed by atomic interactions 
(especially the isotropically produced protons). The probability of even a 100 MeV proton to undergo a nuclear 
reaction is only several percent. The low energy neutrons on the other hand have no effective atomic interactions 
and propagate until a nuclear interaction occurs or they escape from the media. Yet, the low energy protons and 
other light ions produced in tissue through neutron collisions contribute significantly to biological injury and must 
be represented [Foelsche et al. 1974]. This was the basis of the studies of Heinbockel and co-workers [2003] starting 
with Clowdsley et al. [2000, 2001] using the NUCFRG2 nuclear database [Wilson et al. 1987a, 1995, 1998]. 

Development of an algorithm for solving equation (1) or (3) is expedited by approximating the double 
differential cross sections as 

v m (E,E',Sl,n') = a nn j(E,E')6(il - SI') + a nnh (E, E ')6(Sl + SI'). ( 9 ) 

Unlike the straight-ahead approximation that reduces equation (3) to a Volterra equation, this bi-directional equation 
(9) results in a new approximate Fredholm equation that can be solved by a perturbative expansion or other means as 
examined by Heinbockel and co workers. Whereas the straight-ahead solution accounts for leakage only at the distal 
surface, this bi-directional solution provides estimates of the albedo neutron leakage. Although improvements in 
quality of solution can be judged by direct comparison of the differences in the bi-directional solution and the 
straight-ahead values, the overall quality as judged by other codes are limited by the use of differing nuclear 
databases. Still, among the MC codes available, it has been shown that the bi-directional code agrees with MC codes 
to the extent that MC codes agree among themselves in slab geometry [Heinbockel et al. 2011a,b, Slaba et al. 2011, 
2013]. Considered in the subsequent section is the transport in finite objects in which lateral leakage can play an 
important role. 

3D Marching Procedures 

A sequence of approximations that spans the formalism from the simple straight-ahead approximation to 
more complex propagation algorithms is now developed [Wilson et al. 2014a,b,c, 2015]. The governing equation in 
the CSDA ignoring Coulomb scattering is given by equations (1) and (3). The double differential cross sections are 
found to be closely approximated by a forward component and an isotropic component given by 

* jk (E,E'MSl') = <t m JE,E') / 4t r + cr jkJor {E,E\Sl^) , (10) 

where the first term is isotropic and associated with lower-energy particles produced including target fragments 
defined by 
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"Kiso&E') = 2 Jl B * ik (E,E',n,n')dn\ 


(ii) 


where 2nB represents the backward hemisphere and the second term in equation (10) is highly peaked in the forward 
direction and is associated mainly with direct quasi-elastic events and projectile fragmentation products [Wilson 
1977, Wilson et al. 1988b]. The isotropic cross section is given by equation (11). The corresponding forward 
component is found from 


= a jk (E,E',n,n') - a jktiso (E,E') / 47 r . ( 12 ) 

With the demonstrated success of the bi-directional approximation, the solution is likewise divided into forward and 
isotropic components where the isotropic neutron term is treated as a perturbation. The first order term is given by 
the following equations 


B l^fJxME)} = ^ra jkJor (E,E')t kJor (xME')dE' 


(13) 


Whereas the full straight-ahead cross section of equation (8) appears in the ion transport of equation (13), only the 
forward propagating components appear in the transport of neutrons in the equation given by 


a jkJor 


( E,E ') 


f^a jk (E,E',n,n')<m' tj * n 

f 4 }jkjor(E,E',n,n')dn' , i = n 


(14) 


Note that Q enters into equation (13) as a parameter and can be evaluated according to the boundary 
condition that we set as the incident fluence from direction Qo at the boundary as shown in Fig. 2. The coupling to 
the isotropically produced neutron field, treated as a perturbation, is given by 


[n* V + a n (E)}d> n .JxME) = f 4 ^ nn (E,E\n,n')t n .Jx^\E')dn'dE^ 

+E S"^iso(^E\n^ 0 ^(x,n 0 ,E')dE' . 

k 

Once the forward fluence is evaluated for direction Qo over the domain of x and E, one must yet solve for the 
induced isotropic neutron fluence within the shield configuration given by equation (15). Note that the last term in 
equation (15) is an isotropic source of neutrons from collisions of the forward propagating fluence given as 


it 



Fig. 2. Geometry of 3D marching procedure. 
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( 16 ) 




for which the forward fluence (j) k f or (x : El 0 ,E') is a function of the penetration depth z (Fig. 2) along the direction 

Qo. This is denoted by introducing x n i S0 [ z ( x )^^ o>^] = w ^re z(jc) is the penetration depth 

along Qo to the point jc. Equation (15) is rewritten as 


[O . V + a„ (£)R , iso (*, O, £) = f E °° f a nn (E,E',n,Sl ’)</-„ {x,£l\E')dSl' dE' 


' E J 47r 

+X B ,«,[< a! )> n > n o>- B ] 


(17) 


Note that equation (17) also provides a basis for introducing anisotropic source terms as well. Whereas the forward 
propagating solution has its source fixed by the boundary condition, the isotropic solution (having no inbound 
fluence) is driven by internal sources generated by the collisions of the inbound fluence as modified by the forward 
propagator in penetrating to x along direction Qo. In evaluating the penetration depth, one must account for 
differences in the material composition that determines the attenuation in reaching x from the boundary in addition 
to the strength of the neutron source term at location x. 

It is clear that Q enters equation (17) as a parameter and can be solved along an arbitrary number of rays used 
to construct a numerical representation of cj) niso (x,£l,E) whose Q dependence comes from the distance to the 

boundary t( Q) and the penetration depth z along Qo. The isotropic neutron source of equation (17) effectively 
decouples the incident direction Qo from the neutron direction of propagation Q. To solve equation (17), the neutron 
fluence, cj) niso (x,tt,E) , must be known at every x from all directions Q' approaching x from adjacent locations. 

Although the use of approximation (9), based on the assumption that losses along the ray Q are compensated by 
diffusion from adjacent rays, is a demonstrated useful assumption in a semi-infinite slab [Wilson et al. 2014a, b], it is 
an approximate overestimate in any finite object. Still, using equation (9) in evaluation of equation (17) seems a 
reasonable approximation in most human rated vehicles leading to an extension of transport along an arbitrary ray 
Q as a solution to the bi-directional equations 

[«• v + a n (E)]cP n .JxME) = J^°° a nnf (E, E ')<j> niso (x, fl,E')dE' 

+ f„°nnA E , E 'KiJ X -n,E')dE' (18) 

+x n , lso Hx),n,n 0 ,E ] , 

where the backward moving field satisfies a similar equation obtained by replacing Q by -Q in equation (18) as 

[-n*V + a n (E)]X t Jx,-n,E) = f E a nnJ (E,E')ct> n!lS0 (x-n,E')dE' 

+fz*nn, b ( E > E 'KiJxME')dE' (19) 

+X„ )iso Hx),-n,n 0 ,E} . 


These isotropically produced neutrons produce isotropic light ions in collisions with the media resulting in a light 
ion fluence given approximately by 


st.v^x.a,E)-X± Sj(E) 


4> u Jx,n,E) = ( 20 ) 
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where (j) n iso {x,£l, E) is the isotropically produced neutron fluence. Equations (18) and (19) are solved 

simultaneously for no additional external fluence incident on the media directed along Qo (see Fig. 2). Note that the 
implicit assumption of equations (18) and (19) is that lateral diffusion from the ray Q is fully compensated by 
diffusion from adjacent rays. This has been demonstrated as a reasonable assumption for slab geometry [Heinbockel 
et al. 201 la, b; Slaba et al. 2011, 2013; Wilson et al. 2014a, b,c, 2015]. The lateral leakage along Qo in this case is 
found to a first order approximation by solving equations (18) and (19) along each Qi between the ray limits {-£(- 
QO, £(Q0} for a given set of N directions Q,. 


Simple Homogenous Geometry 


Simple geometry and source orientation are chosen to allow efficient MC simulation for comparison with 
the calculation of quantities with the above formalism. In this way, the solution methodology is verified and 
experience is gained in the formalism in preparation for connecting to the general geometric descriptions generated 
by the engineering design process for which MC simulation is computationally difficult. In one recent example, 195 
CPU-days were required for a single evaluation for incident protons in relatively simple geometric models of the 
International Space Station (ISS) [Ersmark et al. 2007]. Even the development of a simplified model at each step in 
the design process would disrupt design optimization if 195 CPU-days for each design option were required [Wilson 
et al. 2004]. Still, the MC methods play an important role in testing the development of highly efficient 
computational methods as has been demonstrated elsewhere [Wilson et al. 1991, 2005, 2014a, b,c, 2015]. 

Consideration is given first to a 40 g/cm 2 homogeneous slab of aluminum with sides extending to infinity. 
A subset of the slab (cube) is also considered, as shown in Fig. 3. The transport within the slab has leakage from the 
cube subset as shown in the upper figure, compensated by lateral leakage into the subset as shown. The Webber 
representation of the February 23, 1956 solar particle event (SPE) [Webber 1966] is assumed incident from above 
(an exponential momentum spectrum with momentum parameter po = 100 MV and 10 9 protons/cm 2 above 30 MeV). 
A convergence test is performed on the slab lateral dimensions with 3DHZETRN for N = 22 rays [Wilson et al. 
2014a,b,c, 2015] and consider the length of the sides required to minimize lateral leakage at the detector locations 
indicated in Fig. 3. Results are shown in Fig. 4 at depths of 35 g/cm 2 and 40 g/cm 2 . It is clear that the impact of 
lateral leakage on the neutron fluence at detector locations is minimal at 500 g/cm 2 and larger. For comparison, 
results from Geant4 version 9.4.6 [Agostinelli et al. 2003] are also shown. Removing the 40 g/cm 2 cube from the 
slab, as shown in the lower half of Fig. 3, leaves the leakage from the cube that is no longer compensated by 
diffusion from the exterior slab into the cube. Note that the neutron collisional source given by equation (16) within 
the cube is the same while embedded within the slab as when exposed singularly. 

In Fig. 4, the 3DHZETRN neutron fluence at detector locations in the boxes have increases at the lowest 
energies with increasing lateral dimensions up to 500+ g/cm 2 . These increases within the present calculation are in 
reasonable agreement with Geant4 above 0.3 MeV as seen in the figure. At the lower energies, where absorptive 
processes are minimal, the leakage along the long transport paths parallel or nearly parallel to the front and back slab 
faces are underestimated in the present formalism, resulting in the increase in fluence at these low energies. 
Improved propagation methods are required to correct this problem. 


^ ^ ^ ^ ^ ^ ^ 
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Fig. 3. Aluminum slab geometry (top) and associated cube geometry (bottom). The slab is 40 g/cm 2 thick with 
infinite lateral dimensions. The cube is 40 g/cm 2 on all sides. The red circles indicate detector locations at which 
quantities are evaluated. The green arrows indicate general directions of leakage. 
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Fig. 4. Convergence test for solution in slab geometry at depths of 35 and 40 g/ cm 2 . 


The calculated fluence using the current method within the 40 g/cm 2 cube at 35 g/cm 2 and 40 g/cm 2 depths 
is shown in Fig. 5 in comparison with three MC codes (Geant4 version 9.4.6 [Agostinelli et al. 2003], FLUKA 
[Fasso et al. 2005, Battistoni et al. 2007] and PHITS [Sato et al. 2006, 2013]. All four codes use differing nuclear 
databases (see Appendix A of Wilson et al. [2014a,b]), which accounts for some of the differences in Fig. 5. In 
addition, the current version of 3DHZETRN assumes the cross sections consist of only straight-ahead and isotropic 
components in the present formalism for the neutron source distributions. The differences among the MC codes is 
assumed mainly due to differing nuclear databases while much of the difference of 3DHZETRN lies in both the 
nuclear database and the forward/isotropic assumption. Herein, attempts are made to separate the effects of the 
forward/isotropic assumption from that of the nuclear database limitations. 



This is accomplished, in part, by looking at the neutron leakage ratio given by 


r(x,E) = 1 — 


object 

<l>n,slab( X ’ E ) ’ 


( 21 ) 
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where <fi n oh j ect (x, E) is the neutron flux obtained at a depth in a finite object, and is the neutron flux 

obtained at the corresponding depth in the semi-infinite slab geometry. The leakage ratio denotes the fractional 
leakage from the location v for neutrons of energy E through near lateral boundaries. The leakage defined by 
equation (21) from the cube is shown in Fig. 6 for 3DHZETRN in comparison with results from three MC codes 
(Geant4, PHITS, and FLUKA). Whereas the fluence estimates of the four codes can differ significantly as seen in 
Fig. 5, the leakage ratio varies only by ten to fifteen percent among the codes as seen in Fig. 6. Also note that the 
vertical axis of Fig. 6 (and subsequent plots showing leakage factors) is in linear scale instead of log scale. Using the 
leakage ratio removes some of the dependence on the magnitude of the neutron production and brings into clearer 
focus the transport processes of scattering and absorption. 



Fig. 6. Comparison of MC evaluated leakage from a cube at 35 and 40 g/cm 2 depths. 


The high leakage indicated in the 3DHZETRN solution below 0.3 MeV results from the neutron fluence 
overestimate along rays nearly parallel to the front and back slab faces in slab geometry as shown in Fig. 4. It is 
clearly seen that this "apparent" error in diffusive losses at low energies of the present code results in an 
overestimate of leakage from the cube and requires improved representation of diffusive losses through the faces in 
slab geometry as was discussed above in connection to Fig. 4. The adequacy of the forward representation in the 
straight-ahead approximation is seen at the highest energies for this geometry. At intermediate energies, the errors in 
leakage are mainly associated with the isotropic angular distribution assumed, and improved angular dependence in 
the neutron source (see equation (16)) is required. 

Another geometry of interest is a spherical subset of the slab as shown in Fig. 7. Unlike the cube, the 
neutron collisional source within the sphere by itself differs from the distribution within a sphere embedded in the 
slab. This occurs since portions of the sphere while embedded in the slab are shielded by the exterior portion of the 
slab. 







Fig. 7. Aluminum slab geometry (top) and associated sphere geometry (bottom). The slab is 40 g/cm 2 thick with 
infinite lateral dimensions. The sphere diameter is 40 g/cm 2 . The red circles indicate detector locations at which 
quantities are evaluated. The green arrows indicate general directions of leakage. 
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While leakage effects from the cube are unambiguous, leakage from the sphere is obscured by the 
differences in the neutron source distribution; still, leakage as defined by equation (21) remains a useful concept. 
The neutron source distributions for energies above 1 MeV in the cube and sphere are shown in Fig. 8. The 
calculated neutron fluence at 35 and 40 g/cm 2 from the present method is shown in comparison with that of the three 
MC codes in Fig. 9. It is first noted that the current solution method is in better agreement with the MC codes in 
spherical geometry as seen in comparing Fig. 9 with Fig. 5. 


1956 WebberSPE ^ ^ 


1956 WebteSPE 



x g/cm 


x g/crrf 


Fig. 8. Isotropic neutron source (E > 1 MeV) distribution in units of neutrons/(g-event) induced by the Webber SPE 

spectrum in cube and spherical geometry. 



Fig. 9. Particle spectra at detector locations 35 g/cm 2 (left pane) and 40 g/cm 2 from top of aluminum sphere 
(diameter 40 g/cm 2 ) exposed to Webber SPE spectrum. 


The results of the present calculation for leakage at 35 and 40 g/cm 2 within the sphere are shown in 
comparison with results from the three MC codes in Fig. 10. All of the codes have differences in nuclear models (the 
HZETRN models are discussed in Wilson et al. [1991], and the MC models are discussed in Appendix A of Wilson 
et al. [2014a, b]). Differing nuclear models are the main source of differences observed in Fig. 10. The present 
calculation uses the straight-ahead approximation for the forward components at the highest energies for which no 
lateral leakage occurs and is the primary source of leakage differences between the cube (Fig. 6) and the sphere (Fig. 
10) above 100 MeV [Wilson et al. 2014a,b]. Clearly the straight-ahead approximation is more appropriate along the 
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central ray of the cube and less appropriate in the sphere. The overestimate of leakage below 0.3 MeV follows from 
errors in diffusion through the slab faces along rays parallel (or nearly parallel) to the facial plane. Improved 
treatment of the nuclear models below 10 MeV, where quasi-elastic scattering is negligible, and above 100 MeV, 
where the assumed straight-ahead approximation introduces error, will require formal improvements to the present 
calculations. 
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Fig. 10. Comparison of MC evaluated leakage from a sphere at 35 and 40 g/cm 2 depths. 


To further examine 3D effects, the leakage of the cube and sphere in the present code is compared with MC 
codes that have better representation of low energy scattering and the high-energy angular dependent cross sections 
(relative to the straight-ahead approximation used herein). MC (Geant4, FLUKA, PHITS) and 3DHZETRN derived 
results are shown in Figs. 5, 6, 9-14. Whereas good agreement of the leakage factor for the cube above 100 MeV is 
obtained (see Fig. 6), the failure of the straight-ahead approximation in the sphere is apparent in Fig. 10. Note that 
the leakage of the sphere exceeds that of the cube in the present formalism and that the main differences are above 
100 MeV. The source of this high-energy difference results from the slight angular divergence from the forward 
direction that is of less importance along the cube centerline compared with leakage along the centerline of the 
sphere. 



Fig. 11. Leakage factors for cube and sphere at 35 g/cm 2 and 40 g/cm 2 evaluated by the 3DHZETRN ( N = 22) code. 
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Surprisingly, the leakage factor for the cube and sphere are largely similar in the 3DHZETRN formalism as 
seen in Fig. 11. Hence, the geometric shape in this case has only modest effects in the present formalism within 
similar sized objects. This is an advantage in that the neutron spectrum depends weakly on geometric detail. An 
interesting feature of the present formalism is that the cube exhibits greater leakage than the sphere at intermediate 
energies (1 to 100 MeV). This results from the isotropic source term which has the same distribution in the cube but 
is significantly modified in the sphere as shown in Fig. 8 and places the intense neutron source region closer to the 
detectors and hence appears as reduced leakage as defined by equation (21). 

This weak dependence on geometry is somewhat borne out in the Geant4 results shown in Fig. 12, except 
that the sphere exhibits greater leakage as intuition would suggest. The greater leakage expected for the sphere 
above 100 MeV is likewise given by the Geant4 code. Similar results are found for the PHITS code as seen in 
comparing Fig. 13 with Fig. 12. A surprising result of the present benchmark is the results from the FFUKA code 
seen in Fig. 14. Although the excess leakage in the sphere above 100 MeV is clearly due to the straight-ahead 
approximation, the intermediate energy excess leakage is due to the differences within the FFUKA nuclear model. 
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It is clear from the above results that improvements in the transport model in which the first order 
scattering angular dependence is treated as an isotropic source lies in the following: the low energy transport 
requires higher order corrections to scattering, and the isotropic source contributions need to be expanded to include 
more complex angular factors extending up to higher energies than currently represented. Further complications are 
added in solving the Boltzmann equation in inhomogeneous objects. The next section examines inhomogeneity in a 
tissue sphere shielded by aluminum. 

Simple Inhomogeneous Geometry 

In this section, simple inhomogeneous slab geometry is considered by substituting a 30 cm (30 g/ cm 2 ) thick 
layer of ICRU tissue [ICRU 1993] within 7.41 cm (20 g/cm 2 ) of aluminum shielding above and below. The cube 
subset becomes a box of 44.82 cm on each side. Both geometries are shown in Fig. 15. The alternate spherical 
geometry will consist of an ICRU tissue equivalent sphere (30 cm diameter) surrounded by an aluminum shell with 
thickness 7.41 cm (20 g/cm 2 ). As in the case of the homogeneous configurations, convergence of the neutron 
fluence in the slab is considered first; with exposure of this arrangement from above by the Webber SPE. Results are 
shown in Fig. 16. 
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Fig. 15. Slab geometry (top) and associated cube geometry (bottom). The slab is 30 g/cm 2 (30 cm) of ICRU tissue 
with 20 g/cm 2 (7.41 cm) of aluminum above and below and infinite lateral dimensions. The cube has the same 
thickness of ICRU tissue and aluminum and has lateral dimensions of 44.82 cm. The red circles indicate detector 
locations at which quantities are evaluated. The green arrows indicate general directions of leakage. 
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Fig. 16. Convergence test in an inhomogeneous slab composed of a 30 cm ICRU tissue slab shielded above and 
below by 7.41 cm of aluminum exposed to the Webber SPE spectrum. Results shown at 0 cm (left pane) and 30 cm 
(right pane), corresponding to the top and bottom detector locations, as shown in Fig. 15, respectively. 


The E >1 MeV isotropic neutron source for the rectangular box and shielded sphere (discussed later in this 
report) are shown in Fig. 17. One sees clearly that the cube neutron source is distributed along the z-axis as one 
would expect for the slab configuration. The abrupt change in crossing from the aluminum shield into tissue is 
clearly seen near the top interface. A similar transition occurs at the top of the tissue sphere as well. 

The calculated neutron fluence generated in the inhomogeneous cube derived from the 3DHZETRN and 
MC codes is shown in Fig. 18. Note that the 3DHZETRN tends to have good overall agreement with the MC 
(Geant4, PHITS, FLUKA) codes at least to the extent the MC codes agree among themselves. The neutron leakage 
as seen from the detector locations on both the top and bottom of the tissue phantom is shown in Fig. 19 as evaluated 
by 3DHZETRN and MC codes. Unlike leakage from the uniform aluminum cube, the tissue phantom of the 
inhomogeneous cube provides significant absorption so that neutron leakage is limited. In this case, there is 
reasonable agreement between 3DHZETRN and the MC simulations. 
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Fig. 17. Isotropic neutron source (E > 1 MeV) distribution in units of neutrons/(g-event) induced by the Webber SPE 
in the rectangular box and shielded sphere. The intensity (color) scale has been adjusted to more clearly show the 

shadow cast by the tissue sphere. 
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Fig. 18. Particle spectra at detector locations 0 cm (left pane) and 30 cm (right pane) from top of the tissue box 
(thickness 30 cm) shielded by aluminum above and below (thickness 7.41 cm) exposed to the Webber SPE 

spectrum. 



Fig. 19. Leakage from an inhomogeneous cube of aluminum/tissue/aluminum as seen at the center of the top 

interface (x = 0) and the lower interface (x = 30). 


Consider now a 30 cm diameter ICRU tissue sphere shielded by a 7.41 cm spherical shell of aluminum, as 
shown in Fig. 20. The change in calculation due to the inhomogeneity is that the penetration depth must evaluate the 
aluminum depth separate from the depth in the tissue sphere, and the solution along any direction Q, will include 
propagation thru both media from the outer limits of the object along the ray ±Q, to the evaluation point. In this 
example, the solution for neutron fluence is evaluated at the top and bottom of the ICRU sphere that serves in this 
work as proxy for the astronaut (ICRU phantom, see ICRU [1993] and Wilson et al. [2014c, 2015]). The neutron 
fluence is evaluated in the tissue and aluminum interface at the top and bottom of the tissue sphere using both 
3DHZETRN and MC methods (Geant4, PHITS, and FLUKA). In order to control statistical fluctuations, the MC 
simulations require a relatively large volume over which fluences are scored. This implicitly integrates physical 
variations found locally near the detector volumes. To provide a consistent comparison, the 3DHZETRN fluence is 
obtained as an average of point results over the same volume (a sphere of 0.25 cm radius). The neutron fluence at 
top and bottom of the tissue sphere are shown in Fig. 21 where reasonable agreement with MC results are obtained. 
The corresponding neutron leakage ratio is shown in Fig. 22. 
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Fig. 20. Slab geometry (top) and associated sphere geometry (bottom). The slab is 30 g/cm 2 (30 cm) of ICRU tissue 
with 20 g/cm 2 (7.41 cm) of aluminum above and below and infinite lateral dimensions. The sphere is an ICRU tissue 
sphere with diameter 30 g/cm 2 (30 cm) surrounded by 20 g/cm 2 (7.41 cm) aluminum. The red circles indicate 
detector locations at which quantities are evaluated. The green arrows indicate general directions of leakage. 




Fig. 21. Particle spectra at detector locations 0 cm (left pane) and 30 cm (right pane) from top of tissue sphere 
(radius 15 cm) surrounded by an aluminum shell (thickness 7.41 cm) exposed from above by the Webber SPE 

spectrum. 


Conclusions 

The focus of the present study was neutron production and propagation. It was shown that the neutron 
differential fluence solution above 1 MeV in slab geometry converges well to match the Geant4 simulation in 
homogeneous materials. At lower energies, an error is introduced by the bi-directional approximation ( N = 2) that 
assumes no additional diffusive losses for the propagation of the isotropic neutron component. The losses through 
the slab faces in transport along rays that are parallel or nearly parallel to the slab faces are underestimated (assumed 
zero). As a result, the neutron leakage ratio for the sphere and box is overestimated at low energies. A correction for 
losses along these long paths near lateral boundaries is required. The convergence for solution in the slab geometry 
is greatly improved when a highly absorptive tissue layer is introduced as an inhomogeneous layer. Above 100 
MeV, the 3DHZETRN generally underestimates the leakage due the straight-ahead assumption and is easily seen in 
the neutron leakage ratio of the sphere where angular divergence is expected to play an important role. On the 
centerline of the cube, the straight-ahead approximation is more correct. These errors will be partly corrected by 
treating the angular dependence in greater detail. 
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Fig. 22. Leakage from an inhomogeneous sphere of aluminum/tissue at the top interface (left pane) and lower 

interface (right pane). 


The codes for space design are still incomplete in both the physical description and the corresponding 
implementation of the computational procedures. Even so, the current development is state-of-the-art in 
computational procedures for spacecraft design. The present version of 3DHZETRN has some inherent limitations 
(aside from the simplified geometry) that carry over from the forward/backward transport methods in which the 
main share of light ions are still treated in the straight-ahead approximation as coupled to the forward propagating 
neutrons. Only the "assumed" isotropically produced neutrons are first transported in 3D without the coupling to the 
isotropically produced light ion components. Once the isotropically produced neutrons are evaluated, the final 
coupling to the light ions is implemented (without further nuclear reactions) in a full 3D manner. Hence, only a 
fraction of the light ions have a 3D treatment resulting in transport effects of unquantified magnitude. Future work 
will require extension to complex geometries and treatment of the cross-coupling of 3D neutron/light ion transport. 
Nonetheless, the present expanded treatment of 3D effects show great promise in the path to a fully 3D transport 
algorithm. 
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